Ultracold bosons in lattices with binary disorder 
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Quantum phases of ultracold bosons with repulsive interactions in lattices in the presence of quenched dis- 
order are investigated. The disorder is assumed to be caused by the interaction of the bosons with impurity 
atoms having a large effective mass. The system is described by the Bose-Hubbard Hamiltonian with random 
on-site energies which have a discrete binary probability distribution. The phase diagram at zero temperature 
is calculated using several methods like a strong-coupling expansion, an exact numerical diagonalization, and a 
Bose-Fermi mapping valid in the hard-core limit. It is shown that the Mott-insulator phase exists for any strength 
of disorder in contrast to the case of continuous probability distribution. We find that the compressibility of the 
Bose glass phase varies in a wide range and can be extremely low. Furthermore, we evaluate experimentally 
accessible quantities like the momentum distribution, the static and dynamic structure factors, and the density 
of excited states. The influence of finite temperature is discussed as well. 

PACS numbers: 03.75.Lm, 03.75.Hh, 67.85.Hj 



C/3 



I 

o 
o 



> 
m 

O 
oo 
O 



I. INTRODUCTION 

The remarkable experimental control over ultracold atomic 
gases in optical lattices acquired in recent years U 12. HI B HI 
I2] has opened up completely new lines of investigation in the 
field of strongly correlated quantum systems. One of these are 
quantum phase transitions (QPTs) of ultracold atoms in opti- 
cal lattices. These fascinating phenomena are caused by the 
interplay of quantum tunneling, atomic interaction and dis- 
order. In contrast to other condensed-matter systems where 
quantum phase transitions can also take place, optical lattices 
provide a unique possibility to control the disorder which can 
be created by several methods. Truly random potentials with 
continuous disorder distribution can be created using laser 
speckles iSiSlii which leads to random contributions to 
the tunneling amplitudes as well as on-site energy shifts. In 
addition to that, the atomic interaction energies can be made 
random ifTHl . if the lattice loaded by cold atoms is placed near 
a wire inducing a spatially random magnetic field II12I1 . Disor- 
der with discrete probability distribution can be created intro- 
ducinga second atomic species strongly localized on random 
sites fisi [Til [Tsll which leads only to random shifts of the on- 
site energies. With the aid of incommensurate lattices one can 
make the tunneling amplitudes and the on-site energies quasi- 
random 

Until recently, studies of QPTs in cold atoms were deal- 
ing with continuous disorder distributions of different types. 
QPTs in the presence of disorder with discrete probability dis- 
tribution were studied only for interacting electrons in binary 
alloys fl^. The role of this type of disorder in QPTs of cold 
bosons starts to become also a subject of research {20I1 . In the 
present work, we shall investigate QPTs of ultracold bosons 
in a lattice with disorder which is created by the interaction 
with impurity atoms localized at random lattice sites. The un- 
derlying theoretical model is the Bose-Hubbard Hamiltonian 
with random on-site energies according to a binary probability 
distribution. The problem was recently addressed in the work 
by Mering and Fleischhauer ll20ll . who employed the DMRG- 
method to study the system. Our analysis will be based either 



on an exact numerical diagonalization for sufficiently small 
systems or a Bose-Fermi mapping for the case of hard-core 
bosons in a one-dimensional lattice. These exact methods will 
be applied to obtain the phase diagram at zero temperature. 
Moreover, the role of finite temperature will be studied. In ad- 
dition, we evaluate experimentally accessible quantities, such 
as the momentum distribution, the static and dynamic struc- 
ture factors and the density of excited states. 

The paper is organized as follows. In Sec. [Ill we describe 
the theoretical model of the Bose-Hubbard Hamiltonian with 
binary disorder, studied in the present work. In Sec. [nil the 
physics of the Mott-insulator (MI) phases is explained and 
their phase-boundaries are calculated by employing perturba- 
tion theory with respect to the hopping amplitude. SectionlTVl 
deals with exact numerical calculations of the many-particle 
ground states for small lattices. In Sec.|V] we consider a one- 
dimensional system in the limit of strong interaction, which 
is exactly solvable through the Bose-Fermi mapping. A sum- 
mary of the work is given in Sec.lVIl 



II. HAMILTONIAN 

We consider a system of ultracold interacting bosons in 
a d-dimensional hypercubic lattice described by the Bose- 
Hubbard Hamiltonian 



(iJ> 



(1) 



where J is the tunneling matrix element for the nearest lattice 
sites, U is the on-site atom-atom interaction energy, and fi 
is the chemical potential. Throughout the paper, we will be 
dealing with repulsive interaction, i.e., U > 0. We assume 
periodic boundary conditions. The annihilation and creation 
operators, ai and Aj, obey the bosonic commutation relations. 
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The disorder described by the random terms ei is assumed 
to be created by the presence of impurity atoms located at 
fixed random positions (quenched disorder). This type of dis- 
order is diagonal in the sense that it does not lead to any contri- 
butions to the hopping term. If it has a "fermionic" character, 
i.e., if there is at most one impurity at each lattice site, the 
probability distribution of on-site energies ei is given by 

pie)^PoS{e) + {l-po)S{e-U') , po = {L-L')/L , (2) 

with < Po < 1, where U' is the boson-impurity interaction 
energy. L is the number of lattice sites and L' is the number 
of impurities, po = 1 or U' = corresponds to the pure case. 
The binary disorder distribution (|2]i implies that the system 
under consideration remains invariant under the transforma- 
tion Po ^ 1 — Po, U' — [/'. This can be easily seen if we 
subtract U' /2 from the on-site energies e,, i.e., if we make the 
replacement e — > e + U' /2 on the r.h.s. of Eq. (|2|. Therefore, 
it is enough to consider the case U' > 0. 

III. STRONG-COUPLING EXPANSION 

It is known that the Bose-Hubbard model with continu- 
ous disorder distribution possesses a rich phenomenology of 
phases resulting in a nontrivial phase diagram ll2lll . For the 
present case of binary disorder, we determine the zero temper- 
ature phase diagram next. The physics of the MI phases can 
be understood and their phase boundaries be calculated with 
a good accuracy by treating the hopping term in the Hamilto- 
nian ([l]i as a perturbation f23\. This can be done in arbitrary 
spatial dimensions and, in this section, we will not impose any 
restrictions with respect to the dimensionality d. For the bi- 
nary disorder distribution (|2]i it is convenient to consider the 
entire lattice as consisting of two disconnected sublattices Cq 
and Ci. The sublattice Ci consists of the L' potential wells 
which are shifted by ei — U' with respect to the wells of 
the other sublattice Cq. In what follows, it is assumed that 
< Po < 1- One may then define the local chemical poten- 
tials of the sublattices Co and Ci as /i and p — [/', respectively. 
The special case of a pure lattice can be retrieved in the limit 
U' -^0. 

In the limit of vanishing hopping (J = 0), the states of the 
system are characterized by the occupation numbers of each 
lattice site. In the ground state corresponding to the MI, we 
have no bosons at each site of the sublattice Co, and ui bosons 
at each site of the sublattice Ci, where no and ni are the 
smallest non-negative integers larger than or equal to and 
(/i — U')/U, respectively. In general, ni < no and the total 
number of bosons N is given by N ~ no{L — L') + niL' . We 
are interested in the thermodynamic limit N oo, L oo, 
L' oo, where N/L as well as po remain finite. The par- 
ticle and hole excitations in the limit of infinite lattices are 
localized within the Lifshitz rare regions ll25ll which consist 
of infinitely large connected regions of either sublattice Co or 
Ci depending on the disorder parameter U' and the chemi- 
cal potential /i. This allows us to work out the boundaries of 
the MI regions by generalizingthe method of strong-coupling 
expansion developed in Ref . |26[| . Due to the infinite extent 
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FIG. 1: The boundaries of the MI phases obtained by the strong- 
coupling expansion in d = 1 (dotted), 2 (dashed), 3 (solid) for 
U/U' = 0.7 (top) and for U/U' = 1 (bottom). 

of the Lifshitz rare regions, the boundaries of the MI regions 
should not depend on po- This dependence appears as a finite 
size effect which we do not consider in this section. 

The region fi < ~2dJ corresponds to a vanishing particle 
number no = '^i = 0. The intervals 7io — 1 < n/U < no, 
7io = 1, . . . , [U' /U], with [. . . ] denoting the integer part, cor- 
respond to the MI with ?ii = 0. The lowest-energy particle- 
hole excitations are created by transferring one atom among 
the lattice sites of the sublattice Co- The energy gap for the 
creation of this excitation equals U at J = 0. The MI phases 
are enclosed in the interval i^i^ino) < fi < /.Jp(no), where 

/ip(n) = Un-2dJ{n+l) (3) 

+ —n[di5n + A)-4d^{n + l)] 

J3 

+ T7:^n{n + l)\-8d^{2n+l) + d^{25n+U) 
-4d(2n+l)] + 0(j^) , 

^^(ji) ^ U{n - 1) + 2dJn 

J2 

- —in+l)[d{5n+l)~-M'^n] 

+ —n{n + l)\8d^{2n+l)^d\25n + ll) 

+Ad{2n+1)] + 0{J^) , 

are the boundaries of the MI phase in the pure case ll26ll . 

The next intervals no — 1 < fx/U < no, where no = 
\U' /U] + 1, . . . , are split into two subintervals. In the lower 
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subinterval no — 1 < ijl/U < ng — 1 + {U' /U}, where {. . . } 
denotes the fractional part, ni ~ no — [U' /U] — 1. The lowest- 
energy particle-hole excitation of this state can be created by 
transferring one atom from the sublattice Cq to the sublat- 
tice Ci. As a result of this transfer, the energy of the initial 
state is increased by U {U' /U}. Perturbative calculations in 
the thermodynamic limit show that it is located in the interval 
Hh{no) < fJ. < U' + fJ.p{ni) of the phase diagram (Fig.[T]). 

In the upper subinterval no — 1 + {W /U} < ji/U < no, 
ni = no — [U'/U]. The lowest-energy particle-hole excitation 
of this state can be created by transferring one atom from the 
sublattice Ci to the sublattice Cq. As a result of this transfer, 
the energy of the initial state is increased hy U ~ U {U' /U}. 
According to the perturbation theory, it is located in the inter- 
val U' + iih{ni) < ji < /ip(?7o), where iih{n) and fJ,p{n) are 
given by Eq. (O. 

The two subintervals exist only if {[/'/[/} does not vanish, 
otherwise the lower subinterval disappears. The upper one 
then becomes extended from fi — {no ~ 1)U to fi ~ noU. 
In this case, it is energetically more favorable to create the 
particle-hole excitations by transferring one atom among the 
lattice sites of the sublattice Co and the MI phase is enclosed 
within the interval /i;i(no) < fi < iip{no). 

MI phases with equal occupation numbers no = ni = n 
exist only for < U' /U < l,i.e., [U'/U] ^0 and {U'/U} = 
U' /U . They are located within the intervals no — 1 + U'/U < 
ijl/U < no. 

The MI regions for U'/U = 0.7 and for U'/U = 1 are 
shown in Fig.[T] In the case U'/U = 0.7, we have only split 
MI regions with no = 1,2, ... , where the lower and upper 
parts correspond to ?ii ~ no — I and ni ~ no, respectively. 
In the case U'/U = 1, there are only nonsplit MI regions with 
n-i = no — 1. With the increase of the dimensionality d, the 
MI regions become smaller. 

Since we are dealing with a disordered system, one can ex- 
pect the existence of the Bose-glass (BG) phase ||2l|). How- 
ever, the method of strong-coupling expansion in its present 
form does not allow us to detect the corresponding regions 
on the phase diagram. It does not give the opportunity to in- 
vestigate the temperature-dependent effects either. Therefore, 
other methods are needed in order to study the complete phase 
diagram of the system. They are employed in the next sec- 
tions, where we consider only one-dimensional lattices. 



IV. EXACT DIAGONALIZATION 

In this section we study zero-temperature properties of the 
system by means of exact numerical diagonalization of the 
Bose-Hubbard Hamiltonian along the lines of Ref. IztIi . For 
this we determine first the boundaries of the regions in the 
{fi, J) plane corresponding to different total particle numbers 
N. This requires calculations of the ground-state energies Ej^ 
of the Hamiltonian ([T]i for different N and can be done exactly 
with the aid of iterative numerical solvers for sparse matrices 
of large dimensions. 

The results of these calculations for a small one- 
dimensional lattice are presented in Fig. |2l The solid lines 
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FIG. 2: The regions with occupation numbers A'^ — 0, . . . , 14 
obtained by the exact diagonalization of the Hamiltonian (TJ for 
Z/ = 10, I/' = 6 for one disorder realization characterized by the 
Fock state |0011101011> for U'/U = 0.7 (a) and for U'/U = 1 
(b). 



indicate the boundaries fi^ = En — E^-i of the regions 
with different occupation numbers N of the lattice. These cal- 
culations are performed for i = 10 and L' ~ 6 and for one 
disorder realization, where the spatial distribution of impuri- 
ties is described by the Fock state 1 00 11101011). The maxi- 
mal number of bosons is 15. The lowest line is the boundary 
between = and = 1, the next one is the boundary 
between = 1 and N ^ 2 and so on. 

In the case U'/U = 0.7 (Fig. |2^), the regions with the 
occupation numbers N~L~L' — A, N~L = 10, 
N = 2L — L' = \A appear to be larger than the others indi- 
cating that there are MI phases for these occupations. In spite 
of the large contribution of the finite-size effects, the shape 
of these regions is in good agreement with the results of the 
strong-coupling expansion, see Fig.[Tl 

As it was discussed in the previous section, the MI phases 
with integer filling factors disappear \fU'>U and only the 
MI phases with incommensurate fillings remain. This behav- 
ior can be also seen in Fig.|2j5. 

With the increase of the number of lattice sites L, the 
boundaries of the regions with different occupation number 
are closer to each other and in the thermodynamic limit they 
should densely cover the whole J) plane except the MI 
regions. In order to really see this transition to the thermody- 
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namic limit as well as to determine the boundaries of the MI 
regions, one has to vary the number L of lattice sites and the 
number A'^ of bosons in a wide range which is difficult here 
because the dimension of the bosonic Hilbert space grows ex- 
ponentially with N and L. 

The fraction of the total atom number N which are 
in the superfluid phase (superfluid fraction) can be calculated 
with the aid of Peierls phase factors. They have to be intro- 
duced in the Hamiltonian ([T]) by means of the replacement 
a\dij^i — !■ al^fli+ie**^. By calculating the free energy F/v (0) 
for some small value of (fs, the superfluid fraction is deter- 
mined as imin 



0^0 



fjv ((/.)- Fat (0) 



(4) 



The limit — > can be calculated exactly if the complete 
solution of the eigenvalue problem is known (see Ref. iIztIi 
and the discussion in Sec. IV Al ). However, all the numerical 
solvers for large sparse matrices allow efficient calculations 
only of a small number of the eigenstates. This is the rea- 
son why the superfluid fraction defined by Eq. Q is usually 
worked out for some small but nonvanishing value of (p. 

The behavior of for different N is shown in Fig. [3] for 
the same fixed disorder realization as above. In these calcula- 
tions, (f) = 0.01. In the case U' /U = 0.7, the superfluid frac- 
tion vanishes not only for fillings which allow the existence of 
the MI phases (N = 4, 10, 14) but also for all the others if the 
hopping parameter J is small enough. This suggests a phase 
transition from the superfluid into the BG phase. The case 
U' /U = 1 looks different. The superfluid fraction vanishes 
only for N = 1,2, 3, 4, 14 and remains finite for all the oth- 
ers, even for small values of J/U. This suggests that the BG 
phase exists in the extended regions of the (//, J) plane corre- 
sponding to low fillings but is strongly suppressed for higher 
fillings. 

Since exact diagonalization can be performed only for very 
small lattices, it is difficult to provide a satisfactory descrip- 
tion of the phase diagram of the system. Finite-size effects 
can be much better controlled in the case of hard-core bosons, 
which is treated in the next section. 



V. BOSE-FERMI MAPPING 

We consider the hard-core limit of infinitely strong repul- 
sion U ^ cx) in a one-dimensional lattice, which is exactly 
solvable via the Jordan- Wigner transformation ll29[[30ll 



ai = exp I itt'^CjCj | q 



(5) 



where c; and cj are the fermionic annihilation and creation op- 
erators. Under this transformation the Hamiltonian ([ij takes 
the form 



H = -J 



E(4 



Ci+l 



^Y^{e,~^i)clc.,. (6) 





FIG. 3: Superfluid fraction vs A'^ and J obtained by exact di- 
agonalization for L = 10, L' = 6 for U'/U = 0.7 (a) and for 
U' /U — 1 (b). The points for different values of A'^ are connected in 
order to guide eyes. 



Periodic boundary conditions for bosons are equivalent to the 
requirement 



CL+i = exp 




(7) 



which implies periodic boundary conditions for fermions if 
the number of particles is odd, otherwise one should use the 
corresponding antiperiodic boundary conditions in the Hamil- 
tonian 

The A^-particle eigenstates of the Hamiltonian ^ can be 
constructed from the L single-particle eigenstates as 



(8) 



with the eigenenergies e^, a = 1, . . . , L. Here, T is the trans- 
lation operator. Its action on the one-particle Fock state 
|10 . . . 0) results in the shift of the particle's position by i lat- 
tice sites. This allows to treat much larger lattices compared to 
the case of soft-core bosons considered in Sec.|IV] However, 
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one still cannot avoid numerics because the analytical solu- 
tion of the single-particle eigenvalue problem in the presence 
of disorder is not known. 

In the infinite-?/ limit, the occupation numbers of the indi- 
vidual lattice sites can be or 1 implying that the maximal 
number of atoms N cannot be larger than L. The state with 
N = Lis always a MI as no hopping can take place any more 
and a non-trivial treatment of the MI phase in the hard-core 
limit is possible only for N < L. 



A. Bose-Einstein condensation and superfluidity 

Before starting the discussion of the remaining part of the 
phase diagram, some remarks on the Bose-Einstein condensa- 
tion (BEC) and superfluidity of hard-core bosons in ID are in 
order In the absence of disorder, the spatial correlation func- 



tion in the limit \i ~ j\ ^ oo decays as \i — 

The presence of disorder makes this decay faster, i.e., there is 

no off-diagonal long-range order and BEC llsill . 

The absence of BEC does not exclude in general the su- 
perfluidity. As it was shown in Ref. |[32ll for the Gaussian 
disorder, the system of one-dimensional soft-core bosons is 
in the delocalized (superfluid) state, if the correlation func- 
tion {a\aj) for large distances decays slower than \i — 
which leads to the divergence of the localization length. Since 
the correlation function of hard-core bosons with disorder de- 
cays faster than |i — the criterion of Ref. lf32ll excludes 
the existence of the superfluid phase. The question is whether 
this remains true for the binary disorder. 

The superfluid fraction is defined by Eq. The limit cf) 
can be calculated making use of the perturbative approach 
of Ref. iIztIi . For hard-core bosons in ID, this leads to the 
following general expression at T = 0: 



L N 



i=l a=l 
L N 

- ± y y^— 

N ^ ^ Ea-en 
Q=Af+l/3=l " 

L 2 



(9) 



1=1 

where (^^(L + 1) = (— In the clean case, the 
second term in Eq. ^ vanishes and in the thermodynamic 
Umit we get = sisi", where n = N/L fH. This result 
does not depend on J. 

The results of numerical calculations for the binary disor- 
der are shown in Fig.|4]for J/U' = 1 for two different lattice 
sizes. Here and later on the overline indicates the disorder- 
averaged quantity. For a small lattice size L = 100, fs ap- 
pears to be finite. However, when increasing Llo L = 200, fs 
approaches zero. For smaller values of J /U' , smaller lattice 
sizes L are sufficient in order to see that fs indeed vanishes 
in the thermodynamic limit, i.e., we expect only insulating 
phases in our system. For comparison, the corresponding re- 
sult for the clean case is also shown in Fig.|4l see dashed line. 
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FIG. 4: Superfluid fraction for J/U' = 1, L' = L/2, L = 100 
(circles), 200 (solid line) averaged over 400 disorder realizations. 
The dashed line shows the superfluid fraction in the thermodynamic 
limit without disorder, i.e., fs 



sin(7riV/L) 
■kN/L 



see text. 
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FIG. 5: The regions with occupation numbers A'^ — 0,...,L 
obtained by the exact diagonalization of the Hamiltonian l|6} for 
L = 200, L' — 100. The results are averaged over 400 disorder 
realizations. Dashed lines are analytical results for the boundaries of 
the MI regions in the thermodynamic limit. 



Since the hard-core Umit works for U' /U < 1 and J/U 
1, this qualitatively agrees with the behavior of ff for N < 
i = 10 and relatively small J/U calculated by means of exact 
diagonalization in the case of soft-core bosons (Fig.[3h). 



B. Quantum phases 

The boundaries of the Ml-regions at zero temperature can 
be easily calculated analytically in the thermodynamic limit 
following the general treatment of Sec.|III]for d = 1, U' < U. 
The region /i < —2 J contains no atoms. The MI phases with 
non-integer filling factors are located in the interval /i/j (uq) < 
fi < U' + /ip(ni), where uq ~ 1, rii ~ 0. In the limit 
U ^ oo, all the terms in Eq. ^ proportional to Jp with p > 
1 vanish and we get explicitly 2J < /.i < U' — 2 J, where 
J < [/'/4. U' + fih{ni) < ^ < /Up (no) with uq = ni = 1, 
i.e., U' + 2J < jj, < oo corresponds to the MI with N = L. 
These results are shown in Fig.|5]by dashed lines. 

The same boundaries as well as the boundaries of the re- 
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gions with different occupation numbers /i n can be calculated 
numerically by means of exact diagonalization. are equal 
to the single-particle eigenenergies en because the ground- 
state energy of N non-interacting fermions is the sum of N 
lowest single-particle eigenenergies. The region /i < ei con- 
tains no atoms. The MI phase with L — L' atoms is located 
in the interval sl-L' < yU < el-L'+i- If A' > ^l, we have 
the MI with N = L. The boundaries ei, . . . , are shown 
by solid lines in Fig. |5]for L = 200. The lines ei, El-L', 
£l are outside of the corresponding regions deter- 
mined in the thermodynamic limit, similar to what is seen in 
Fig. 121 They come closer to the results obtained in the limit 
L ^ oo if the lattice size L is increased. 

The distribution of lines /ijv = £n{J), N ~ 1, . . . , L, in 
Fig.|5]is very inhomogeneous which leads to the fact that the 
compressibility of the system varies in a rather wide range. 
This characteristic feature remains preserved for larger lattices 
as well and it is easier to see it in the behavior of N{fj,) which 
is given by 



L 

E 

a=l 



fiSa), /(£) 



1 



exp [(e - + 1 



(10) 

The plots N{^) at T = are shown in Fig.|6h. The central 
plateaus (A^ ~ L — L') around n/U' = 0.5 which exist for 
J/U' < 0.25 correspond to the MI phase. Quasi-plateaus of 
N{ii) which exist at any values of J/U' correspond to the re- 
gions on the {J, /i)-diagram with low density of lines fiN{J). 

In order to clarify the physical interpretation of different 
parts of the phase diagram, we have calculated the time- 
dependent Green's function G(r) = (ai(T)aJ(O)) defined 
for r > which determines the density of states of the 
single-particle excitations as well as the superfluid suscepti- 
bility [21J. The plots of G(r) as a function of the imaginary 
time r are shown in Fig.]?] In the MI phase, G{t) at large t is 
an exponential function of r. This also holds for N = L — L', 
J/U' < 0.25. If the number of particles remains the same but 
J/U' is increased, the exponential decay of G(t) is replaced 
by the 1 / r-law indicating that now we are in the BG phase 
(Fig.|2]main) according to Ref. ll2lll . 1/r-law is also observed 
for other particle numbers corresponding to the quasi-plateaus 
of N{^) (inset of Fig. |7]) which allows to interpret them as 
belonging to the BG phase, in spite of the fact that the com- 
pressibility is extremely small. 

The density of states for the single-particle excitations can 
be determined in terms of the Fourier transformed single- 
particle Green's function G(i!;) as p{E) = -^lmG{E) B . 
For hard-core bosons, it takes the form 



piE) 



f{Ep)]6iE~Ep+E^) 



(11) 



where Ef3 — Ea are the energies of single-particle excitations 
caused by the transfer of one particle from the energy level Ea 
to the energy level e^. 

Numerical calculations were performed assuming that the 
energy levels have a finite lifetime. The 5-function is approxi- 
mated by a Gaussian with the standard deviation O.OIC/'. The 
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FIG. 6: (Color online) N{p) in the hard-core limit for ksT /U' = 
(a), 0.1 (b) and for one disorder realization which is the same in (a) 
and (b). L = 5000, L' = 2500, J/U' = 0.1 (black), 0.2 (red), 0.3 
(blue). 



disorder-averaged energy dependences at T = are shown 
in Fig. [8] for different values of J. 'p{E) has always a multi- 
maxima structure which stems from the inhomogeneous dis- 
tribution of the single-particle eigenenergies (see Figs.|5]|6^). 

The density of excited states al E = as a function of 
J is shown in Fig. |9] for = L — L' . p{Q) vanishes for 
J/U' < 0.25, and is different from zero otherwise. For 
other particle numbers, p(0) does not vanish for any J. This 
is directly related to the asymptotic properties of the time- 
dependent Green's function G(t) discussed above. p(0) > 
is a clear-cut signature of the BG ll2ll l24ll . 

At finite temperature, N{p) defined in Eq. ( fTOl i becomes a 
smooth function (Fig.|6}3). The MI plateaus which are clearly 
seen for J/U' = 0.1,0.2 at T = (Fig. |6ji) are smeared 
out and the quasi-plateaus disappear even at rather small T 
indicating that the compressibility never vanishes. Neverthe- 
less, one can say that the MI still exists as long as the gap for 
the creation of particle-hole excitations, which is U' — 4 J in 
the case of hard-core bosons in ID, is larger than the thermal 
energy ksT. This suggests that the upper and lower bound- 
aries of the Ml-region are given by U' — 2 J — kBT/2 and 
2J + kBT/2, respectively, and the crossover line for the MI 
on the (J, T) -diagram is given by ksT^" = U' ~ 4 J. 

If temperature is increased, quantum statistics will become 
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FIG. 7: (Color online) Linear-log plot of G(r) for L = 400 averaged 
over 5000 disorder realizations for = 200 and J /U' — 0.1 (black 
solid), J/U' = 0.2 (red circles) and J/U' = 0.3 (blue stars). Inset: 
same in a log-log plot for J/U' = 0.3 for = 132 (solid) and N = 
200 (dashed) which corresponds to fi/U' ~ and 0.5, respectively. 
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FIG. 8: (Color online) Density of excited states at T = averaged 
over 200 disorder realizations. L = 400, L' = 200, TV = 200, 
J/U' = 0.1 (black), 0.2 (red), 0.3 (green), 0.4 (blue). 



less important and one can expect that the MI phase as well as 
the BG phase is destroyed. Then, a crossover into the normal- 
gas state is expected. Since the thermodynamic properties of 
hard-core bosons are equivalent to those of the ideal Fermi 
gas, the crossover temperature coincides with the Fermi tem- 
perature modified by the disorder. 



C. Experimentally measurable quantities 

The Bose-Fermi mapping allows very detailed investiga- 
tions of different physical properties of the system which can 
be directly measured in exp eriments. We consider first the 
momentum distribution 13411 



exp [ika{l - I')] (ajai,) , 

1,1' 
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0.1 



0.2 0.3 

J/U' 



0.4 0.5 



FIG. 9: Density of excited states at T = and E = averaged over 
200 disorder reahzations. L = 400, L' = 200, N = 200. 



where k is the wavenumber, a is the lattice constant, and 
W{k) is the Fourier transform of the Wannier function W{x) 
for the lowest Bloch band of the lattice potential. The matrix 
elements (a|a;/) of the L x L one-particle density matrix can 
be worked out for Z > I' as a determinant of the (Z— /') x (/—/') 
ToepHtz matrix G''^''' iSEll as 



{ajai') 



-yl-l'-l 



The matrix elements of G'' ' are given by 



G 



{1,1 



(12) 



(13) 



The expectation values (clcj) can be calculated using the so- 
lution of the single-particle eigenvalue problem as 

L 

(ctc,) = V^:(z)^„(j)/(£j. (14) 



The momentum distributions obtained by numerical calcu- 
lations for N = L — L' sae shown in Fig. [lOl In general, 
P{k) = {k)4'{k)) /\W{k)\'^ is an even and periodic func- 
tion of ka with the period 27r. It takes maximal (minimal) 
values at ka = 7rm, m = 0, ±2, ±4, ... (m = ±1, ±3, . . . ). 
With the decrease of the hopping parameter J the spatial cor- 
relations of bosons become weaker which leads to the broad- 
ening of the momentum distribution. As discussed above, for 
N = L — L' the system undergoes a phase transition from the 
BG to MI, where the spatial correlation functions obey power 
and exponential laws, respectively. Therefore, the dependence 
of the momentum distribution on J is expected to be weaker in 
the MI than in the BG, which is demonstrated in Fig.[lO] The 
transition point is seen as a kink in the J-dependence of the 
momentum distribution function at fc = (inset of Fig. [TOb . 
For N ~ {L ~ L')/2, which always corresponds to the BG 
phase, the dependence is almost linear without any kink. 

At finite temperature, the qualitative form of the momentum 
distribution remains unchanged. As shown in Fig.[TT] P(0) is 
a decreasing function of T for large enough J/U' correspond- 
ing to the BG region in the phase diagram at T = 0. The 
largest choice J — OAU' is in the BG-region rather far from 
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FIG. 10: (Color online) Momentum distribution in the hard-core 
limit at r = averaged over 200 disorder realizations for L — 200, 
L' = 100, A'^ — 100. Inset: Maximum of the momentum distribu- 
tion. 
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FIG. 11: (Color online) Temperature dependence of -P(O) for L = 
200, L' = 100, A'^ = 100 averaged over 250 disorder realizations. 
Inset: Log-log plot of P{k) for the same parameters. 



the BG-MI transition point implying that the averaged mo- 
mentum distribution has a tendency to become broader due to 
the influence of the thermal fluctuations. For smaller values of 
J /U' , i.e., closer to the BG-MI transition, -P(O) first increases 
for small T and then decreases further. For the smallest value 
J = O.IL/', the system is deep in the MI phase at T = 0, and 
the intermediate maximum in P(0) as a function of T van- 
ishes again. 

Useful information about the state of the many-body system 
can be obtained with the aid of Bragg spectroscopy il35[[36ll . 
The response of the system to this kind of measurement is 
described by the dynamical structure factor, which is defined 

as 



dt{Ap{k,0)Ap{-k,t))cxp{-iujt) , (15) 



where Ap{k) is the spatial Fourier transform of the density- 
fluctuation operator. In the case of deep lattices it takes the 




FIG. 12: (Color online) Dynamical structure factor in the hard-core 
limit at r = averaged over 600 disorder realizations. Parameters 
are L = 200, L' = 100, = 100, ka = 7r/3. Moreover, J/U' = 
0.1 (black), 0.2 (red), 0.3 (green), 0.4 (blue). 



form 

Ap{k) = loik) ^ (^ajai - {ajai}^ exp{ikal) , (16) 
I 

where 

dx cxp{ikx) \ W{x)\'^ 



Io{k) 



(17) 



with W{x) being the Wannier function for the lowest Bloch 
band. The dynamical structure factor for hard-core bosons 
can be expressed in terms of the single-particle eigenmodes 
as il 



Sik,E) = h\ioikfY. 

a, 13 



X fiea)[l-f{£p)]S{E-ep + e^) . (18) 

This formula resembles Eq. (fTTT ) for the density of excited 
states but has a more complicated structure due to the explicit 
dependence on the mode-functions <^a(0- 

The resulting dependence of the disorder-averaged dynam- 
ical structure factor S{k,E) on energy is shown in Fig. [12] 
for ka = 7r/3. The behavior of S{k,E) is completely dif- 
ferent compared to the case of homogeneous lattices studied 
in Ref. 1381] . It vanishes in the finite interval of E near zero, 
provided that J/U' < 0.25, due to the energy gap in the ex- 
citation spectrum of MI. With the increase of J/U' the gap 
decreases. It disappears completely if J/U' > 0.25 due to 
the transition into the BG phase. S{k, E) is broader in the 
BG phase (J/U' = 0.3, 0.4) than in the MI phase {J/U' = 
0.1, 0.2). Its multi-peak structure is qualitatively related to the 
density of excited states, which is shown in Fig. [8] However, 
the detailed form of the energy dependence of S{k,E), which 
remains preserved for larger lattices as well, is different from 
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p{E) due to the nontrivial contributions of the eigenfunctions 
(paii) in Eq- (flSl ). The same features are observed for other 
values of ka. 

Finally, we consider the static structure factor defined as 

/oo 
S{k,E)dE. (19) 
-oo 

In the case of hard-core bosons, it is given by Eq. ( fTSb . where 
the (5-function is formally replaced by 1. Its J-dependence is 
shown in Fig.[T3]for ka = 7r/3,27r/3 and ka = n. In gen- 
eral, it grows monotonously with increasing J. Interestingly 
enough, we find a kink which corresponds to the MI-BG tran- 
sition point, similar to the behavior of P{Q) in Fig.fTOl 

VI. CONCLUSIONS 

In the present work, we have studied quantum phase transi- 
tions of ultracold bosons with repulsive interaction in a lattice 
with binary disorder. The system is described by the Bose- 
Hubbard model with random on-site energies which follow a 
binary probability distribution. The particular form of disor- 
der is physically realized, e.g., when two species of alkali- 
metal atoms with different masses are loaded in an optical lat- 
tice. The latter is created by counter-propagating laser beams 
which are strongly detuned from the atomic resonance. If one 
of the species has a mass, say, 4 times bigger (as it is realized 
for the combination of ^^Rb and -^'^Na atoms), the tunneling 
of the heavier atoms in a deep enough lattice is suppressed 
by more than 3 orders of magnitude compared to that of the 
lighter ones. Thus, they effectively form immobile impuri- 
ties interacting with the lighter atom species. Larger mass 
differences, like for the combination of ^°K-^Li or ^^Rb-^Li 
and ^■^■^Cs-^Li, render the difference in the tunneling rates 
even more drastic. The impurity atoms then induce an ef- 
fectively quenched disorder potential for the lighter bosonic 
atoms. When their number is less than the number of the 
lattice sites, the probability to find more than one impurity 
at a lattice site can be extremely low. This is garanteed by 




J/U' 

FIG. 13: (Color online) Static structure factor as a function of J in 
the hard-core limit at T = averaged over 200 disorder realizations. 
L = 200, L' = 100, N = 100. 



repulsive interactions between the bosonic impurities and by 
Pauli's exclusion in the case of the fermionic ones. Since the 
impurities are assumed immobile, their statistics does not play 
any role. In fact, the interaction parameters U and U' can be 
tuned over a wider range by the additional use of Feshbach 
resonances. 

To calculate the boundaries of the MI phases in the phase 
diagram for arbitrary lattice dimension d at zero temperature, 
we have applied the method of strong-coupling expansion. We 
have shown that the MI phase exists also for incommensurate 
bosonic fillings, and not only for commensurate ones. Fur- 
thermore, the binary disorder generates additional Mott lobes 
in the phase diagram. 

For ID lattices, we have investigated the superfluidity of 
soft-core bosons at T = for the binary disorder. Due 
to the exponential growth of the dimension of the bosonic 
Hilbert space for increasing boson numbers and numbers of 
lattice sites, exact numerical diagonalization is possible only 
for small lattices, which does not allow to control finite-size 
effects. However, the obtained results for small lattices are in 
good agreement with perturbative results obtained in Sec. Hill 
as well as with the exact results in the hard-core limit, see 
Sec.lVAl 

In the limit of infinitely strong repulsion (hard-core 
bosons), we have performed rather detailed exact studies by 
applying the Jordan- Wigner transformation. This allows to 
considerably reduce the computational complexity since all 
the properties of the strongly interacting system can be deter- 
mined in terms of the solution of the single-particle problem. 
The remaining disorder average can straightforwardly be per- 
formed by standard numerical means. We have shown that 
the binary disorder destroys the superfluidity in the thermo- 
dynamic limit in a similar manner as for Gaussian disorder 
studied earlier. However, in contrast to the case of Gaussian 
or uniform disorder, we have found that the compressibility 
of the EG phase can be extremely low. Several experimen- 
tally measurable quantities such as the momentum distribu- 
tion, the dynamical and static structure factors and the density 
of excited states have been worked out. The MI-BG transi- 
tion can be identified via rather sharp kinks in the functional 
dependence of the maximum of the momentum distribution 
on the tunneling ,/ (Fig.[Tol). Similar kinks occur in the static 
structure factor (Fig. \T3[ . These kinks allow to identify the 
MI-BG quantum phase transition. The energy-dependence of 
the dynamical structure factor yields complementary informa- 
tion about the gap in the excitation spectrum in both phases 
(Fig. \12\ . Given the wide availability of elaborated experi- 
mental techniques, we hope that these predicted features will 
be found in real physical systems of ultracold atoms in the 
near future. 

In the present work, we did not consider the effects of har- 
monic confinement which is normally present in most experi- 
ments with ultracold atoms in optical lattices. One can expect 
that our results will remain valid for shallow traps. On the 
other hand, coexistence of different phases in different spa- 
tial regions can come into play as in the case without disor- 
der lUlQ] depending on the range of values of the local chemi- 
cal potential in the region occupied by the atoms. For instance. 
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in the hard-core limit, there might be coexistence of the BG 
and MI with one boson per site, if J/U' > 0.25. In the op- 
posite case J/U' < 0.25, MI with non-commensurate filling 
can coexist with the previous two phases. Detailed studies of 
the coexisting phases and experimental signatures of the cor- 
responding QPT is a separate problem which requires further 
investigations. 



Acknowledgment 



This work was supported by the SFB/TR 12 of the German 
Research Foundation (DFG). 



[1] W. Zwerger, Advances in Solid State Physics 44, 277 (2004). 
[2] D. Jaksch and P. ZoUer, Ann. Phys. (NY) 315, 52 (2005). 
[3] I. Bloch, Nature Physics 1, 23 (2005). 

[4] O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78 179 (2006). 

[5] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, 
A. Sen(De), and U. Sen, Advances in physics 56, 243 (2007). 

[6] I. Bloch, J. Dalibard, W. Zwerger, Rev. Mod. Phys. (to be pub- 
lished), arXiv:0704.3011 

[7] J. E. Lye, L. Fallani, M. Modugno, D. S. Wiersma, and M. In- 
guscio, Phys. Rev. Lett. 95, 070401 (2005). 

[8] C. Fort, L. Fallani, V. Guarrera, J. E. Lye, M. Modugno, 
D. S. Wiersma, and M. Inguscio, Phys. Rev. Lett. 95, 170410 

(2005) . 

[9] D. Clement, A. F. Varon, M. Hugbart, J. A. Retter, P. Bouyer, 
L. Sanchez-Palencia, D. M. Gangardt, G. V. Shlyapnikov, and 
A. Aspect, Phys. Rev. Lett. 95, 170409 (2005). 

[10] T. Schulte, S. Drenkelforth, J. Kruse, W. Ertmer, J. Arlt, 
K. Sacha, J. Zakrzewski, and M. Lewenstein, Phys. Rev. Lett. 
95, 170411 (2005). 

[11] H. Gimperlein, S. Wessel, J. Schmiedmayer, and L. Santos, 
Phys. Rev. Lett. 95, 170401 (2005); Appl. Phys. B 82, 217 

(2006) . 

[12] P Kruger, Phys. Rev. A 76, 063621 (2007); S. Wildermuth, 
S. Hofferberth, 1. Lesanovsky, E. Haller, L. M. Andersson, 
S. Groth, 1. Bar- Joseph, P. Kruger, J. Schmiedmayer, Nature 
(London) 435, 440 (2005). 

[13] P Vignolo, Z. Akdeniz, and M. P Tosi, J. Phys. B 36, 4535 
(2003). 

[14] U. Gavish and Y. Castin, Phys. Rev. Lett. 95, 020401 (2005). 
[15] B. Horstmann, J. 1. Cirac, and T. Roscilde, Phys. Rev. A 76, 
043625 (2007). 

[16] B. Damski, J. Zakrzewski, L. Santos, P. Zoller, and M. Lewen- 
stein, Phys. Rev. Lett. 91, 080403 (2003). 

[17] L. Sanchez-Palencia and L. Santos, Phys. Rev. A 72, 053607 
(2005). 

[18] L. Fallani, J. E. Lye, V. Guarrera, C. Fort, and M. Inguscio, 
Phys. Rev. Lett. 98, 130404 (2007). 

[19] K. Byczuk, M. Ulmke, and D. Vollhardt, Phys. Rev. Lett. 90, 
196403 (2003); K. Byczuk, W. Hofstetter, and D. Vollhai'dt, 
Phys. Rev. B 69, 045112 (2004); K. Byczuk and M. Ulmke, 
Eur. Phys. J. B 45, 449 (2005); K. W. Kim, J. S. Lee, T. W. Noh, 
S. R. Lee, and K. Char, Phys. Rev. B 71, 125104 (2005); 



N. Paris, A. Baldwin, and R. T. Scalettar, Phys. Rev. B 75, 
165113 (2007). 

[20] A. Mering and M. Fleischhauer, Phys. Rev. A 77, 023601 
(2008). 

[21] M. P A. Fisher, P B. Weichman, G. Grinstein, and D. S. Fisher, 

Phys. Rev. B 40, 546 (1989). 
[22] Similar treatment can be found in Ref. (20ll where the analysis 

was restricted by rf = 1 and not all the possible values of the 

disorder parameter U' were explicitly considered. 
[23] A. A. Abrikosov, L. P. Gor'kov and 1. Ye. Dzyaloshinskii, 

Quantum Field Theoretical Methods in Statistical Physics 

(Pergamon Press, New York, 1965). 
[24] K. V. Krutitsky, A. Pelster, and R. Graham, New J. Phys. 8, 187 

(2006). 

[25] 1. M. Lifshitz, Sov. Phys. Usp. 7, 549 (1965) [Usp. Fiz. Nauk 

83, 617 (1964)]; F Cyrot-Lackmann, J. Phys. C 5, 300 (1972). 
[26] J. K. Freericks and H. Monien, Europhys. Lett. 26, 545 (1994); 

Phys. Rev. B 53, 2691 (1996). 
[27] R. Roth and K. Burnett, Phys. Rev. A 67, 031602(R) (2003); 

ibid. 68, 023604 (2003). 
[28] M. E. Fisher, M. N. Barber, and D. Jasnow, Phys. Rev. A 8, 

1 1 1 1 (1973); W. Krauth, Phys. Rev. B 44, 9772 (1991). 
[29] P Jordan and E. Wigner, Z. Phys. 47, 631 (1928). 
[30] E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. (NY) 16, 407 

(1961). 

[31] A. De Martino, M. Thorwart, R. Egger, and R. Graham, 

Phys. Rev. Lett. 94, 060402 (2005). 
[32] T. Giamarchi and H. J. Schulz, Europhys. Lett. 3, 1287 (1987); 

Phys. Rev. B 37, 325 (1988). 
[33] In the case of half-filling, this reproduces a well-known result 

fs = 2/tt. See, e.g., S. Rapsch, U. Schollwock and W. Zwerger, 

Europhys. Lett. 46, 559 (1999). 
[34] V. A. Kashumikov, N. V. Prokof'ev, and B. V. Svistunov, 

Phys. Rev. A 66, 031601(R) (2002). 
[35] T. Stoferle, H. Moritz, C. Schori, M. Kohl, and T. Esslinger, 

Phys. Rev. Lett. 92, 130403 (2004). 
[36] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation 

(Clarendon Press, Oxford, 2003) and references therein. 
[37] P Vignolo and A. Minguzzi, J. Phys. B 34, 4653 (2001). 
[38] G. Pupillo, A. M. Rey, and G. G. Batrouni, Phys. Rev. A 74, 

013601 (2006). 



